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TWO-DEGREE-OF-FREEDOM FLUTTER SOLUTION FOR A PERSONAL COMPUTER 


David L. Turnock 


SUMMARY 

A computer program has been written for a personal computer which can be 
used to make two-degree-of-freedom (bending and torsion) flutter calculations 
by utilizing two-dimensional Theodorsen aerodynamics. The proqram may be used 
to approximately account for Mach number (compressibility) effects and aspect 
ratio (finite span) effects. This report contains the equations of motion, a 
program listing, user instructions, and test cases. 


INTRODUCTION 

Flutter is a self-excited dynamic oscillation produced by a coupling of 
inertial and elastic forces with aerodynamic forces resulting from elastic 
deformations of an aircraft wing. These oscillations could result in 
significant structural damaoe. Consequently, flutter must be taken into 
account during aircraft design; that is, the wing must be designed so that it 
will not flutter within the operating envelope of the airplane. 

The types of vibrations that an airplane wing could experience in flight 
are illustrated in Figure 1. 

If the wing is disturbed at flight conditions below the flutter boundary, 
the ensuing motion is a damped sinusoidal decay inq oscillation as shown in 
Figure la. This is a stable condition. At the flutter boundary, any 
disturbance will cause the wing to oscillate at a constant amplitude as shown 
in Figure lb. This is a neutrally stable condition. At conditions above the 
flutter boundary, a disturbance will produce a divergent oscillation as shown 
in Figure lc. This is an unstable condition. 

The purpose of this paper is to describe a computer program written for a 
personal, or home, computer that can be used to analyze aircraft wings for 
two-degree-of-freedom (bending and torsion) flutter. The equations of 
Theodorsen and Garrick (references 1 and 2) are implemented in this program. 
The analysis is a two-dimensional, two-dearee-of-freedom, representative- 
section method. Although this method is not new, it is sufficiently accurate 
for many present day aDplications. 

The method employed is adequate for moderate-to-hioh aspect ratio wings 
and unswept or slightly swept wings. However, the analysis has certain 
limitations. The aerodynamics are only applicable at subsonic velocities. 
The analysis as implemented in the computer program does not include the 
effects of concentrated masses such as wing mounted engines and fuel tanks. 
It also does not include the effects of aerodynamic control surfaces such as 
ailerons. These limitations are very important to note. Furthermore, the 
analysis used by the program is only an approximation. Therefore, the results 
should only be used as a guide, not as the final authority to determine 
whether or not a new or old wing design is safe from flutter. 



This report presents documentation of the computer program. Also 
presented are the derivation of the equations, user instructions, program 
listinq, and test cases which serve as illustrative examples. 
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SYMBOLS 

non-dimensional distance from midchord to elastic axis, positive aft 
semi chord, c/2, ft 
chord, ft 

span efficiency factor, function of taper ratio 
fundamental frequency in plunge or bending, Hz 
fundamental frequency in pitch or torsion, Hz 
artificial damping 

structural damping in plunge and pitch, respectively 
vertical deflection of the wing in plunge, ft 

plunge velocity of deflecting wing, ft/sec 

plunge acceleration of deflecting wing, ft/sec^ 

/FT, imaginary coefficient 

reduced frequency, bto/V 

mass per unit length of span, slugs/ft 

dynamic pressure, pv2/2, lb/ft^ 

radius of gyration of the winq, non-dimensional 

full winq span, ft 

distance of center of gravity from leadinq edqe of winq, ft 

distance of elastic axis from leading edge of wing, ft 

distance of center of gravity location aft of the elastic axis, 

non-dimensional 

area of the wing, ft^ 

aspect ratio 

flutter determinant coefficients 
Theodorsen circulation function 
center of gravity location, percent chord 
lift curve slope, per radian 

flutter determinant coefficients 

elastic axis location, percent chord 

real part of Theodorsen's circulation function 

imaginary part of Theodorsen's circulation function 

altitude, ft 

pitch inertia per unit length of span, sluq-ft^/ft 

lift per unit length of span, lb/ft 

lift due to plunge, non-dimensional 

lift due to pitch, non-dimensional 

Mach number 

moment per unit length of span, ft-lb/ft 
moment due to plunge, non-dimensional 
moment due to pitch, non-dimensional 
wing imbalance, sluo-ft. 
velocity, ft/sec 
speed of sound, ft/sec 
equivalent airspeed, ft/sec 
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a angle of attack, wing deflection about elastic axis in pitch, 

radi ans 

a pitch velocity of wing deflection about elastic axis, radians/sec 

a pitch acceleration of wing deflection about elastic axis, 

radians/sec^ 

p Mach number correction factor 

u mass ratio, m/irpb^, nondimensional 

p air density, slugs/ft^ 

w coupled frequency, rad/sec 

u>h fundamental plunge (bending) frequency, rad/sec 

co a fundamental pitch (torsion) frequency, rad/sec 

a variable of substitution, complex frequency 

8r real part of complex frequency 

8 1 imaginary part of complex frequency 

EQUATIONS OF MOTION 

The flutter analysis used in this proqram is a two-dimensional, two- 
degree-of-freedom, representative-section method. The term "two-dimensional" 
applies to the aerodynamics that are utilized in the program. A two- 
dimensional wing has a constant lift and moment alonq the span. "Two- 
degree-of-freedom" means that the wing is free to deform or vibrate in two 
different motions. In this case, the motions are bending and torsion. 
Bending, or plunge, is the spanwise bending of the wing about the root. 
Torsion, or pitch, is the chordwise twisting of the wing about the rotation 
point or elastic axis of a beam-like wing. These motions are assumed to be 
simple harmonic or sinusoidal. "Representative section" means the winq 

characteristics can be represented by the properties and motions of a unit 
section of the wing. This section is chosen to be located at the 
three-auarter span station. The elasticity of the wing can be represented as 
two springs as shown attached to the wing section in Fiqure 2. If the wing is 
given a displacement, the wing will oscillate in simple harmonic motion 
(sinusoidal oscillations) at frequencies associated with the two degrees of 
freedom. Flutter occurs when these two modes of motion become coupled and 
extract energy from the airstream to produce a self-excited oscillation. 

The governing differential equations for a two degree-of-freedom flutter 
problem are derived in reference 3 and repeated here as equations 1 and 2. 
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These generalized forces are functions of the bending deflection h and 
the torsional twist a. These forces are defined explicitly in reference 3 
and repeated here as equations 3 and 4. 


L' = Trpb 3 <o 2 {L h h/b + [L a - (1/2 + alLja] 


(3) 


M' = W pb 4 (o 2 {[M h - (1/2 + a)L h ]h/b + [M a - (1/2 + a)(L o + M h ) (4) 


+ (1/2 + a) 2 L h ]a} 

where L. , L , M. , M are defined by reference 3 as 
h a* h a 

L h = 1 - 2i(l/k)(F + iG) (5) 

L a = 1/2 - i ( 1/k ) [ 1 + 2(F + iG)] - 2(l/k) 2 (F + IG) (6) 

M h = 1/2 (7) 

M = 3/8 - i (1/k) (8) 

a 


These are termed Theodorsen aerodynamics because of their reliance on the 
Theodorsen circulation function. 


C = F + iG (9) 

The actual function is determined by several Bessel functions and other 
complicated mathematics. However, an approximation is given in reference 4 
that is dependent upon the reduced frequency, k. 


C(k) = 1 


0.165 
0.0455 . 
k— 1 


1 


0.335 
O" . 
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Theodorsen assumed several things in the derivation of the aerodynamics. The 
major assumption was that the wing is two-dimensional. From this assumption 
follows the assumption of infinite span. To further simplify the 
aerodynamics, incompressible flow was assumed. However, both finite-span and 
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compressible-flow corrections can be included to determine their effect on the 
flutter characteristics. These corrections are defined in the following text. 


The assumption of an infinite span is that no wing tip exists and the 
lift distribution is constant along the span of the wing. The lift curve 
slope of a two-dimensional wing is actually greater than the lift curve slope 
for a finite aspect ratio wing at the three auarter span location. For this 
reason the infinite aspect ratio assumption causes conservative results. 
However, for finite-span analysis, the equation for estimating the average 
lift curve slope is derived in reference 5 as 


C L = r (11) 

a L 
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where C[_ is C[_ for infinite aspect ratio. 
^*00 ^ 


hence 


L « = T TTT[ — 

a 

00 

ire^AR 


( 12 ) 


where 

L aoo = original value for L a from the Theodorsen aerodynamics 
Ci = 2u from theoretical analysis 

L OCco 

e^ = .85 for untapered wings 

If the aspect ratio is assumed to he infinity as Theodorsen did, the 
corrected Cl q is equal to Cl^. 
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Furthermore, a reduced lift means a reduced moment. Because the moment 
Is equal to the lift multiplied by a distance, the moment is adjusted by the 
relationship 


M 


a 




(13) 


where the value M^ is the original value calculated by the Theodorsen 
aerodynamics. 


The assumption of incompressible flow means that the effects of Mach 
number are not considered in the calculations. The correction factor for 
compressible flow, which can be used below a Mach number of 0.8, is derived in 
Reference 6 as 


0 = 



(14) 


All Theodorsen aerodynamics are multiplied by this correction factor. 


FLUTTER DETERMINANT 

The solution of the flutter equations is described in this section. 

After the aerodynamic forces in equations 3 and 4 are substituted into 
equations 1 and 2 and simple harmonic motion is assumed, the two simultaneous 
equations that must be solved are given by 
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(20) 


E = ur 2 [l - («o o /m? (l+ig a )] - 1/2 (1/2 + a) +.M a 


- L o (l/2 + a) + L h (l/2 + a)2 

Several methods exist for solving the system of equations 15 and 16. A 
popular method, and the one used here, is from Theodorsen and Garrick 
(references 1 and 2). An artificial damping g and complex frequency 8 are 
introduced where 

8 = (u a /») 2 (l+1g) (21) 

This term is introduced into equations 15 and 16 and the A" and F expressions 
become 


A = y[l-(w h /w a ) 2 (l+ig h )n]+l h 


( 22 ) 


E = yr 2 [l-n(l+ig a )] - l/2(l/2+a) + M o - L o (l/2+a) + L h (l/2+a) 2 (23) 


The B and F expressions remain unchanged. 

For a non-trivial solution to exist, the determinant of the coefficients must 
vanish. That is. 


IT - B C = 0 (24) 

Equation 24 is quadratic in 8. The flutter solution occurs when 8 is 
real that is, when the artificial damping vanishes. Because the aerodynamic 
forces are parametric in reduced frequency k, equation 24 must be solved 
repeatedly to find the flutter solution. The individual solutions are usually 
plotted as shown in Figure 3 and the flutter velocity and frequency are 
obtained from the qraph. The flutter velocity is defined as the velocity 
where the damping value is zero. The flutter frequency is the corresponding 
frequency at that velocity. 
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Individual values are found each time equation 24 is solved for a value 


of reduced frequency. These values of velocity and frequency are obtained 
from the two resulting complex solutions by using equations 25 through 30. 

0 

= Q r + iOj 

(25) 

n R 

= (a> a /u>) 2 

(26) 

°I 

= (u a /w) 2 g 

(27) 

g 

= Oj/O r 

(28) 

0) 

aia 

(29) 



V 

= bai/k 

(30) 


Equations 5 through 8, 12 through 14, and 22 through 30 were used as 
shown in the text in the "computerized" flutter solution. Equations 10 and 17 
through 20 were used in a slightly different form due to the use of complex 
mathematics. A complete listing of the program source code is found in 
Appendix A. 


USER INSTRUCTIONS 

Of course, any computer program is useless without the knowledge of how 
to use it. This section explains the input parameters needed by the program, 
discusses how to obtain specific numerical values, and describes various input 
options. A description of the output is also included. 

Input of Structural Properties 


Structural values may be difficult to obtain in the form required by the 
program. However, reference 7 and the following guidelines should help 
determine the values of certain parameters needed for input to the computer 
program. 

Mass per unit length of span, m . This is mass in slugs (1 slug = 32.2 lbs) of 
a one foot section of the wing. Assuming an untapered straight wing with 
constant mass properties throughout the value can be derived as 


mass of wing (slugs) 
s (ft)' - 
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If the wing is tapered, swept, or has non-constant mass properties, the mass 
is determined for a one-foot-wide representative section centered about the 
three-quarter semi-span. Assuming a linear distribution of mass, the mass 
value is given by 


•_ mass of wing (sluqs) 
m - Z s T F t) " 

Center of gravity location . This is the location in percent chord of the 
point at which the wing balances on a knife edge as seen in Figure 4. 


CG% = (x CG /c)xlOO 


Vling moment of inertia about rotation point, I„ . This is the pitch inertia 
o? a one-foot representative section o7 the wing measured in slua-ft 2 . 
Experimental measurement may be done with a bi filar pendulum, (see reference 
7, pp. 30-31.) 

Elastic axis location . This is defined as the point about which the wing 
twists as illustrated in Figure 5. 


EA% = (xea/c) x 100 


Semi chord length, b . This is one half the chord length of the representative 
section. 

Natural bendinq and torsion frequencies f. ,f . These are fundamental 
h a 

frequencies in Hertz of the wing in bending and torsion. These are functions 
of mass, inertia, and stiffness. (See reference 7, pp. 31-33.) 

Structural damping - bending and torsion, gh,g a * This is defined as any 
internal damping of the structure. Values of .Oi to .03 are typical, however, 
for a solution that is conservative, a value of zero may be used. 

Aspect ratio, AR . This is the ratio of span to chord defined as 

AR = s 2 /A 


For straight, untapered wings, this reduces to 

AR = s/c 

For a fully two-dimensional analysis, an infinite aspect ratio should be 
used. Infinity however, cannot be input to a computer because of numerical 
difficulties. Therefore, a very large (but finite) value must be used. As 
very few wings have an aspect ratio greater than twenty, a value that is 
significantly greater than this (100,000,000 for instance) is recommended. 
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In making calculations to determine if a wing is free from flutter, it is 
prudent to assume that the aspect ratio is infinite because this assumption 
tends to cause conservative results. 

Input for Graphic Output 


The graphic input determines the plot resolution on the screen. A large 
scale may be needed for some solutions. However, that scale may be far too 
large for others. A typical plot appears in Figure 6. Guidelines for these 
input values follow: 

Maximum frequency to be plotted . This is usually the nearest multiple of five 
above the torsion frequency. 

Maximum velocity to be plotted . This should be about twice the expected 
flutter speed. 

Maximum G-value . This is the maximum value of artificial damping that is 
plotted. A value of .05 is sufficient for most cases. 

The above input will create two plots, frequency versus velocity and damping 
versus velocity. A typical plot has points off of the page. However, this is 
only undesirable if the flutter solution crossing point is off of the page. 

Intermediate Values Option 


As has been previously described, a flutter solution consists of solving 
the equations at many different reduced frequency values. Two values of 
velocity, coupled frequency and damping are obtained each time the equations 
are solved. If intermediate values are to be printed, then every value will 
be printed. Appendix B illustrates a solution with all values being printed. 

Air Density Option 


The air density that the wing flies at is very important to the flutter 
characteristics. A wina may flutter at one velocity and altitude yet be very 
stable at the same velocity and a higher altitude. Therefore, air density 
must be input to the problem also. Two options exist within the program for 
this purpose. 

Option 1 - One value for air density is input. This corresponds to flying the 
wing at one altitude. 

Option 2 - Several values for air density. The computer will calculate a 
solution at each air density or altitude. 

If option 2 is selected, then the following values must also be input: 

Number of values - This is the number of times the program must calculate a 
flutter solution using a different value of air density each time. 

Interval between values - This is the increment to be added to the previous 
density value for the next calculation. 
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First Value - This is the starting value for the computer to calculate. 


Example. Three values with an interval of .001 and a starting value of .001 
are as follows: .001, .002, and .003. 

Mach Number Option 


If compressibility effects are desired, a Mach number must be input. The 
input number must be less than one. If no compressibility effects are 
desired, a Mach number of zero should be used. 

Program Output 


The following information is output by the program: 

Structural parameters . The structural input parameters are printed so that 
the parameters can be saved and known for the flutter solution. 

Graphic input parameters . The graphic input parameters are printed in order 
to be able to determine the plot size. 

Preliminary values. These values are printed in order to know the Mach Number 
and air density. These values appear before each solution. 

Intermediate values . These are printed only if the intermediate values option 
is "turned on." The values that are printed are the reduced freguency k and 
the corresponding values of wing frequency, velocities, and damping. These 
intermediate values are used to form data points which are then plotted. 

Flutter values . The flutter values are printed along with an appropriate 
message. Because there are two sets of values (one for bending and one for 
torsion), the flutter conditions (frequency, velocity) are the numbers with 
the corresponding damping or g-value that is close to zero. 

The plots . The intermediate points of frequency and damping are plotted 
versus velocity. The points labeled "1" in Figure 3 are one set of 
intermediate points. As was stated before, two solutions existed for each 
reduced frequency. The points labeled "2" are the solutions of the 
determinant for a second value of reduced frequency. 

The damping curves indicate that one mode will cross the axis. The point 
of crossing is the corresponding velocity at which the wing will flutter. The 
other mode will generally become very stable as airspeed increases. 

Matched Point Solution 

Because the air density is an important factor in flutter calculations, 
an aircraft may flutter near sea level yet be perfectly safe at the same 
velocity when flying at ten thousand feet. Due to this fact, a flutter 
boundary is sometimes defined as Mach number M versus dynamic pressure q. 

A matched point is defined as the density condition where the flutter velocity 
is the same velocity determined from the product of the Mach number and the 
speed of sound V a . To find a matched point solution, multiple air densities 
are used in flutter solutions at the same Mach number. The velocities are 
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plotted against air density and a curve is faired through the points. Each 
value of air density corresponds to a particular altitude in the standard 
atmosphere for which the speed of sound also has a particular value. For the 
given Mach number a curve of velocity versus density is aenerated for the 
standard atmosphere. The point at which the two curves cross is the matched 
point. The corresponding velocity and air density are used to calculate 
dynamic pressure. When plotted as dynamic pressure against Mach number, a 
flutter boundary is defined. A complete matched point solution is oiven in 
Appendix C. 

CONCLUDING REMARKS 

A personal computer flutter solution has been written that employs 
modified two-dimensional Theodorsen aerodynamics for two degree-of-freedom 
lifting surfaces. The solution is only adequate for unswept or slightly swept 
wings with no concentrated masses located on the lifting surface. The 
aerodynamics utilized are only adequate for subsonic velocities. Control 
surface flutter is not addressed by the program. These limitations are 
further stressed by noting that the program is only to be used as a guide and 
not the final authority for determining if a wing is free from flutter. 
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APPENDIX A 


Copy of Computer Program Listing 

The flutter solution was programmed on an Apple IIE computer. Printed and 
plotted solutions were obtained with an Epson printer with Graftrax Interfaced 
by a Microbuffer II interface.* The computer was configured with 64K of 
internal memory and the program used less than 40K. This program was written 
in Applesoft BASIC. Obviously, the program will not run on other small 
computers unless it is modified to account for differences in other 
manufacturers BASIC codes. However, the conversion of the code should be 
relatively straightforward in most instances. 

The listing of the program begins on the next page. 


*Use of trade names does not constitute an official endorsement, either 
expressed or implied, by the National Aeronautics and Space Administration. 
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10 REM FLUTTER PROGRAM — MODULAR SUBROUTINES 
20 PR# 1 

30 PRINT CHR* <9)j"lL" 

s PRINT CHR* (9) ; "80N” 

40 PR# 0 
50 GOSUB 1000 
60 GOSUB 2000 
70 GOSUB 3000 
80 GOSUB 3500 
90 GOSUB 4000 
100 GOSUB 4500 
110 PRINT 

: PRINT 

: INPUT "MACH NUMBER ? " ; MCH 
120 FOR L = 1 TO NR 
130 PR# 1 

s PRINT CHR* (12) 

140 PRINT “MACH NUMBER s "{MCH}" AIR DENSITY s "}RHO(L)j" SLUGS/FT'NS" 

150 PR# O 

1 60 GOSUB 500 

170 GOSUB 6000 

180 FOR Q = 1 TO 25 

190 GOSUB 5000 

200 GOSUB 7000 

210 GOSUB 8000 

220 GOSUB 9000 

230 GOSUB 10000 

240 GOSUB 11000 

250 GOSUB 12000 

260 GOSUB 13000 

270 GOSUB 17000 

280 GOSUB 14000 

290 NEXT Q 

300 GOSUB 15000 

310 PR# 1 

320 PRINT CHR$ (9) ; "G2RDL" 

330 TEXT 
340 PR# 0 

350 FOR TD = 1 TO 500 
s NEXT TD 
360 NEXT L 

370 INPUT "DO YOU WANT TO SEE ANOTHER MACH NUMBER ? ";R$ 

380 IF R$ = "Y" THEN 110 
390 END 

500 REM INITIALIZE GRAPHICS 

510 HGR2 

520 HCOLOR= 2 

530 HPLOT 0,0 TO 0,42 

: HPLOT 0,50 TO 0,191 
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540 FOR X = 0 TO 278 STEP 2 
550 HCOLOR= 2 
560 HPLOT X ,42 

■ HPLOT X , 1 17 
570 HCOLOR= 1 
580 HPLOT CX + 1),42 

s HPLOT (X + 1) ,117 
590 NEXT X 
600 HPLOT 275,44 

s HPLOT 279,44 
s HPLOT 275,45 
: HPLOT 279,45 
: HPLOT 277,47 
610 HPLOT 275,119 

s HPLOT 279,119 
: HPLOT 275, 120 , 
s HPLOT 279,120 
s HPLOT 277,122 
620 HCOLOR= 2 

: HPLOT 276,46 
: HPLOT 278,46 
630 HPLOT 276,121 

s HPLOT 278,121 
640 RETURN 

1000 REM DIMENSION VARIABLES 

1010 DIM OMEGAF <2) ,K(126) ,KI (126) ,LH<2> ,LA<2> , MOMENT (2) ,C0EFF<3,2) ,VEL<2 
) ,0(2,126) , KN (2) , KP (2) 

1020 RETURN 

2000 REM SET CONSTANTS 
2010 PI - 3.14159265 
2020 R1 = 0.165 
sR2 - 1 
:R3 = 0.335 
:R4 = 1 
2030 II = O 

: 12 = - .0455 

s 13 = 0 
s 14 = - .3 

2040 RETURN 

3000 REM GET INPUT DATA 
3010 HOME 

3020 INPUT "MASS PER FOOT SPAN (SLUGS) ? M ; MASS 

3030 INPUT " C-G LOCATION <7. OF CHORD) ? ";CG 

3040 INPUT "WING MOMENT OF INERTIA ABOUT ROTATION POINT (SLUG-FT~2) ? "5 
I ALPHA 

3050 INPUT "ELASTIC AXIS LOCATION (7. OF CHORD) ? " ;EA 
3060 INPUT "SEMI -CHORD LENGTH (FEET) ? " ; SEMCH 
3070 INPUT "NATURAL BENDING FREQUENCY (HZ) ? ";0H 
3080 INPUT "NATURAL TORSION FREQUENCY (HZ) ? " ; OA 

3090 INPUT "STRUCTURAL DAMPING - BENDING (7. CRITICAL DAMPING) 7 ";GH 

3100 INPUT "STRUCTURAL DAMPING - TORSION (7. CRITICAL DAMPING) ? ";GT 

3110 INPUT " FULL SPAN ASPECT RATIO 7 " 5 RTIO 
3120 PRINT 

: PRINT 

3130 HOME * 

3140 REM CHECK PARAMETERS FOR CORRECTNESS 
3150 PR# 1 

3160 PRINT CHR$ (12) 

3170 PRINT "INPUT PARAMETERS : " 

3180 PRINT 
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3190 

3200 

3210 

3220 

3230 

3240 

3250 

3260 

3270 

3280 

3290 

3300 

3310 

3320 

3330 

3340 

3350 

3360 

3370 


3380 

3390 


3400 

3410 

3420 

3430 

3440 

3450 

3500 

3510 

3520 

3530 

3540 

3550 

3560 

3570 

3580 

3590 

3600 

3610 

3620 

3630 

3640 

3650 

3660 

3670 

3680 

3690 


: PRINT "MASS PER FOOT SPAN : 11 ; MASS* H SLUGS " 

PRINT 

PRINT "CG LOCATION : ";CG;» 7. CHORD" 

PRINT 

PRINT "PITCH INERTIA : ",-IALPHA;" SLUG-FT~2 " 

PRINT 

PRINT "ELASTIC AXIS LOCATION s " ; EA; " % CHORD " 

PRINT 

PRINT "SEMI -CHORD s ";SEMCH;" FEET" 

PRINT 

PRINT "NATURAL BENDING FREQUENCY s ";0H;" HZ" 

PRINT 

PRINT "NATURAL TORSION FREQUENCY s ";0A;" HZ" 

PRINT 

PRINT "BENDING DAMPING s " ; GH; " X CRITICAL DAMPING" 

PRINT 

PRINT "TORSION DAMPING s ";GT;" ‘/.CRITICAL DAMPING " 

PRINT 

PRINT "FULL SPAN ASPECT RATIO « "jRTIO 

PRINT 

: PRINT 

s PR# 0 

PRINT 

s INPUT "IS THIS CORRECT (Y/N) ? ";A* 

PRINT 
: PRINT 
s PRINT 

IF A$ < > "Y" THEN 3000 

EA = EA / 100 
GH ■ GH / 100 
GT =» GT / 100 
CG = CG / 100 
RETURN 

REM GET PLOT DATA 
HOME 

INPUT "MAXIMUM VELOCITY TO BE PLOTTED (FEET/SEC) ? ";VMAX 
PRINT 

INPUT "MAXIMUM FREQUENCY TO BE PLOTTED (HZ) ? " ; FMAX 
PRINT 

INPUT "MAXIMUM G-VALUE TO BE PLOTTED ? ";GMAX 

PR# 1 

HOME 

PRINT 

: PRINT "MAXIMUM VELOCITY TO BE PLOTTED : ";VMAX;" FEET/SEC" 
PRINT 

s PRINT "MAXIMUM FREQUENCY TO BE PLOTTED s " ; FMAX; " HZ" 

PRINT 

s PRINT "MAXIMUM G-VALUE TO BE PLOTTED s " ; GMAX 
PRINT 
: PR# 0 
PRINT 

: INPUT " IS THIS CORRECT (Y/N) ? ";A$ 

PRINT " " 

IF A$ < > "Y" THEN 3500 

HOME 

PRINT 

s INPUT "ARE INTERMEDIATE VALUES TO BE PRINTED ? ";NS 
PR# 1 
PRINT 
s PRINT 
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3700 PR# 0 

3710 RETURN 

4000 REM READ K-VALUES 

4010 FOR I - 1 TO 25 

4020 READ K ( I ) 

4030 KHI) = 1 / K ( X ) 

4040 NEXT I 

4050 DATA 10.00,6.00,4.00,3.00,2.00,1.50,1.20,1.00,0.80,0.66,0.60,0.56, 
0.50,0.40,0.30,0.20,0. 16,0.12,0.10,0.08,0.06,0.04,0.025,0.01,0.001 
4060 RETURN 

4500 REM READ IN OR SET RHO VALUES 
4510 HOME 

4520 PRINT "INPUT OPTIONS ON RHO s " 

: PRINT 

4530 PRINT " 1. SINGLE RHO-USER DECIDED" 

4540 PRINT " 2. MULTIPLE RHO'S-USER DECIDED " 

4550 PRINT 

: INPUT "OPTION ? "}CH 
4560 IF CH = 1 THEN 4650 
4570 PRINT 

: INPUT "HOW MANY RHO VALUES DO YOU WANT TO INPUT "jNR 
4580 PRINT 

s INPUT "INTERVAL BETWEEN VALUES ? ";RI 
4590 PRINT 

: INPUT "BEGINNING RHO VALUE ? " ; BR 
4600 DIM RHO (NR) 

4610 FOR L = 1 TO NR 

4620 RHO (L) = BR + (L - 1) * RI 

4630 NEXT L 

4640 RETURN 

4650 NR = 1 

4660 INPUT "RHO VALUE 7 " ; RHD < 1 ) 

4670 RETURN 

5000 REM CALCULATE THEODORSON VALUES 
5010 LI = 12 / K (Q) 

:L2 = 14 / K (Q) 

5020 H = (Rl * R2 + II * LI) / ( (R2 A 2) + (LI ~ 2)) 

5030 J = (R3 * R4 + 13 * L2) / ( (R4 ^ 2) + (L2 ~ 2)) 

5040 N = 1 - H - J 

5050 HI = (R2 * II - Rl * LI) / ( (R2 2) + (LI ~ 2)) 

5060 J I = (R4 * 13 - R3 * L2) / ( (R4 A 2) + (L2 ~ 2)) 

5070 NI = 0 - HI - JI 

5080 LH ( 1 ) = 1 + 2 * NI * KI (Q) 

5090 LH (2) = - 2 * N * KI(Q> 

5100 LA ( 1 ) = 1 / 2 + 2 * KI (□) * NI - 2 * N * KI (Q) 2 

5110 LA (2) = - KI (Q) - 2 * KI (Q) * N - 2 * (KI (Q) ~ 2) * NI 

5120 LR ( 1 ) = LA(1> 

: LR (2) = LA (2) 

5130 FOR Z = 1 TO 2 

5140 LA ( Z) = (1 / (1 + 2 / (.85 * RTIO) ) ) * LR(Z) 

5150 NEXT Z 

5160 MOMENT (1) =3/8 

5170 MOMENT (2) = - KI (Q) 

5180 FOR Z = 1 TO 2 

5190 MOMENT (Z) = LA(Z) / LR(Z) * MOMENT ( Z > 

5200 NEXT Z 

5210 BETA = 1 / SQR (1 - MCH ~ 2) 

5220 FOR X = 1 TO 2 
5230 LH (X) = LH ( X) * BETA 
5240 LA (X) = LA (X ) * BETA 
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5250 MOMENT (X) = MOMENT (X) * BETA 
5260 NEXT X 
5270 RETURN 

6000 REM SET INTERMEDIATE CONSTANTS 
6010 MU = MASS / (PI * RHO(L) * SEMCH * 2) 

6020 IF MU >4 THEN 6050 
6030 PR# 1 

s PRINT "WARNING - LOW MASS RATIO MAY INVALIDATE EQUATIONS " j 
6040 PR# 0 

6050 XALPHA « 2 * CG - 2 * EA 

6060 RALPHA2 = I ALPHA / (MASS * (SEMCH ^ 2) > 

6070 PS I = MU * RALPHA2 
6080 PHI = -1/2+2* EA 

6090 DLTA = (OH / OA) ~ 2 
6100 TAU = MU * XALPHA 
6110 RETURN 

7000 REM SET COEFFICIENTS 

7010 COEFF (1,1) = MU * DLTA * PSI * ( 1 - GH * GT) 

7020 COEFF (1,2) = MU * DLTA * PSI * (GT + GH) 

7030 COEFF (2,1) = MU * DLTA * (PHI / 2 - MOMENT ( 1 ) + LA(1) * PHI - LH(1> * 

PHI ^ 2 - PSI - LA (2) * GH * PHI + GH * LH(2> * PHI ^ 2 + MOMENT (2) * 

GH) + PSI * (LH (2) * GT - LH(1) - MU) 

7040 COEFF (2,2) - MU * DLTA * (PHI * LA(2) - PHI ~ 2 * LH(2) - MOMENT (2) + 

PHI / 2 * GH - MOMENT ( 1 ) * GH + LA(1) * PHI * GH - GH * LH(1) * PHI ~ 

2 - PSI * GH) + PSI * ( - LH(1) * GT - LH(2> - MU * GT) 

7050 COEFF (3,1) = PSI * LH ( 1 ) + MOMENT ( 1 ) * LH ( 1 ) + MU * PSI - MU * PHI / 

2 + MU * MOMENT (1) - MU * LA ( 1 ) * PHI + MU * LH ( 1 ) * PHI * 2 - MOMEN 

T (2) * LH (2) - TAU 2 + 2 * LH ( 1 ) * PHI - 1 / 2 * LA ( 1 ) - TAU / 2 - 

TAU * LA<1) 

7060 COEFF (3,2) = LH<2) * (PSI - PHI / 2 - LA ( 1 ) * PHI + LH(1) * PHI ~ 2) 

+ (LH ( 1 ) + MU) * (MOMENT (2) - LA(2) * PHI + LH(2) * PHI * 2) + LH(2 

) * PHI * (TAU + LA(1) - LH(1) * PHI) - (LA (2) - LH(2) * PHI) * (TAU 

+ 1 / 2 - LH ( 1 ) * PHI) 

7070 RETURN 

8000 REM SET 2ND INTERMEDIATE VALUES 

8010 DR = COEFF (2,1) ~ 2 - C0EFF(2,2) * 2 - 4 * C0EFF(1,1) * C0EFF(3,1) + 

4 * COEFF (1,2) * COEFF (3,2) 

8020 DI = 2 * COEFF (2,1) * C0EFF(2,2) - 4 * COEFF ( 1 , 1 ) * C0EFF(3,2) + 4 * 
COEFF (1,2) * COEFF (3,1) 

8030 SI = SQR ( ( - DR + SQR (DR A 2 + DI 2) ) / 2) 

8040 SR = DI / (2 * SI) 

8050 RETURN 

9000 REM SET 3RD INTERMEDIATE VALUES 

9010 Y — ( ( — COEFF (2,1) + SR) * 2 * COEFF ( 1 , 1 ) + ( - C0EFF(2,2) + SI) * 

2 * COEFF (1,2) ) / (4 * (COEFF (1,1) ~ 2 + C0EFF(1,2) ~ 2)) 

9020 Z — ( ( — COEFF (2,1) - SR) * 2 * C0EFF(1,1) + ( - C0EFF(2,2> - SI) * 

2 * COEFF (1,2) ) / (4 * (C0EFF(1,1) * 2 + C0EFF(1,2) ^ 2)) 

9030 U = (2 * COEFF (1,1) * ( - C0EFF(2,2) + SI) - ( - C0EFF(2,1) + SR) * 

2 * COEFF (1,2)) / (4 * (COEFF (1,1) ~ 2 + COEFF (1,2) 2)) 

9040 W * (2 * COEFF (1,1) * ( - C0EFF<2,2> - SI) - ( - C0EFF(2,1) - SR) * 

2 * COEFF (1,2)) / (4 * (C0EFF<1,1> * 2 + C0EFF(1,2) 2)) 

9050 RETURN 

10000 REM SOLVE FOR FLUTTER FREQUENCIES 
10010 IF Y < 0 THEN 10100 
10020 IF Z < 0 THEN 10100 
10030 OMEGAF(l) « DA / SQR (Y) 

10040 OMEGAF (2) - OA / SQR (Z) 

10050 RETURN 
10100 TEXT 

10110 PRINT "NEGATIVE NUMBER UNDER SQUARE ROOT WHEN TRYING TO SOLVE FOR 
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FLUTTER FREQUENCY "5 
10120 PRINT 

1 PRINT "CALCULATIONS ABORTED " ; 

10130 END 

11000 REM SOLVE FOR FLUTTER VELOCITY 

11010 VEL(l) = SEMCH * OMEGAF ( 1 ) * 2 * PI t K(Q) 

11020 VEL(2) = SEMCH * OMEGAF (2) * 2 * PI / K(Q> 

1 1030 RETURN 

12000 REM SET 4TH INTERMEDIATE VALUES 
12010 Cl = (OA / OMEGAF (1)) /v 2 

12020 C2 = (OA / OMEGAF (2) > ~ 2 

12030 RETURN 

13000 REM SOLVE FOR G-VALUES 
13010 G(1,Q> “ U / Cl 
13020 G(2,Q) = W / C2 
13030 RETURN 

14000 REM PRINT OUT PRELIMINARY DATA 
14010 IF N* < > "Y" THEN 14040 

14020 PR# 1 
14030 PRINT 

14040 PRINT "K-VALUE : " 5 K ( Q ) 

14050 PRINT 

s PRINT "FREQUENCIES (HZ) s "; OMEGAF ( 1 >; " " 5 OMEGAF (2) 

14060 PRINT 

: PRINT "VELOCITIES (F/S) s “jVEL(l>}" ";VEL(2> 

14070 PRINT 

s PRINT "G-VALUE 1 ">G(1,Q)|" ";G(2,Q> 

14080 PRINT 
: PRINT 
14090 RETURN 

15000 REM DWELL ON FLUTTER 

15010 MR = 0 

15020 RL = O 

15030 FOR E = 1 TO 2 

15040 FI = 0 

15050 FOR M = 1 TO 25 

15060 SI = G(E,M> 

:S2 = G(E,M + 1) 

15070 IF SI < O AND S2 > O AND FI = 0 THEN GOSUB 15400 
150B0 NEXT M 
15090 NEXT E 

15100 IF KN(1> > KN (2) THEN R = 1 
15110 IF KN (2) > KN ( 1 ) THEN R = 2 
15120 KL = KN(R) 

: KH = KP (R) 

15130 FOR Q = 26 TO 125 
15140 K (Q) = (KL + KH) / 2 
15150 KI(Q) = 1 / K(Q) 

15160 GOSUB 5000 

15170 GOSUB 7000 

15180 GOSUB 8000 

15190 GOSUB 9000 

15200 GOSUB 10000 

15210 GOSUB 11000 

15220 GOSUB 12000 

15230 GOSUB 13000 

15240 GOSUB 17000 

15250 GOSUB 16000 

15260 IF MR = 1 THEN 15310 

15270 IF G(1,Q> > O AND G(2,Q) > O THEN KH = K(Q) 
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15280 IF G ( 1 , Q) < 0 AND G(2,Q) < 0 THEN KL - K <Q) 

15290 IF (G(1,Q) > 0 AND G<2 f Q) < 0) OR (G(1,Q) < 0 AND 6(2,0) > 0) THEN 
KH = K(Q) 

15300 NEXT Q 
15310 RETURN 
15400 FI = 1 

sKNCE) ■ K (M) 

:KP(E) = K (M + 1) 

: RETURN 

16000 REM PRINT DWELLING DATA 
16010 IF N* < > "Y" THEN 16040 

16020 PR# 1 
16030 PRINT 

16040 PRINT "K-VALUE : " ; K (Q) 

16050 PRINT 

: PRINT "FREQUENCIES i “ ; OMEGAF < 1 ) ; " OMEGAF (2) 

16060 PRINT 

s PRINT "VELOCITIES : " 5 VEL < 1 ) 5 " ";VEL(2) 

16070 PRINT 

: PRINT "G-VALUES : ";G(1,Q)5" *';G(2,Q) 

160B0 IF RL = 1 THEN 16170 

16090 IF ABS (VEL ( 1 ) - VA) < =* 1 OR ABS <VEL<2> - VB) < » 1 THEN 161 

20 

16100 VA = VEL ( 1 ) 

: VB = VEL (2) 

16110 GOTO 16190 
16120 RL = 1 
: MR « 1 
16130 PR# 1 
16140 PRINT 

16145 IF N$ = "Y" THEN PRINT CHRS (12) 

16150 PRINT "THE AIRCRAFT WILL FLUTTER NEAR THE FOLLOWING VALUES " 

16160 GOTO 16030 
16170 PRINT 

16180 PRINT GMAX * ( - 1);" STABLE DAMPING UNSTABLE ";GM 

AX;" O FREQ (HZ) " ; FMAX; 

16190 RETURN 


17000 

REM 

PLOT 

DATA 




17010 

XI 

= 

INT 

(VEL ( 1 ) 

/ 

VMAX * 

280) 

17020 

X2 

= 

INT 

(VEL (2) 

/ 

VMAX * 

280) 

17030 

Y 1 

= 

INT 

(G ( 1 ,Q) 

/ 

GMAX * 

75) 

17040 

Y2 

= 

INT 

(G(2,Q) 

/ 

GMAX * 

75) 

17050 

Z1 

= 

INT 

(OMEGAF ( 1 ) 

/ FMAX 

* 42) 

17060 

Z2 

=5 

INT 

(OMEGAF (2) 

/ FMAX 

* 42) 

17070 

Y 1 

= 

117 ~ 

Y1 





:Y2 = 117 - Y2 
: Z 1 = 42 — Z 1 


s Z2 a 42 - Z2 

17080 IF XI > 279 THEN 17160 
17090 IF Z1 < 0 THEN 17140 

17100 IF XI = (2 * INT (XI / 2)) THEN GOTO 17120 
17110 HCOLOR= 1 

: GOTO 17130 
17120 HCOLOR= 2 
17130 HPLOT X1,Z1 

17140 IF Y1 < 50 OR Y1 > 191 THEN 17160 

17150 HPLOT X1,Y1 

17160 IF X2 > 279 THEN 17240 

17170 IF Z2 < O THEN 17220 

17180 IF X2 = (2 * INT (X2 / 2)) THEN GOTO 17200 
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17190 HC0L0R= 1 

: GOTO 17210 
17200 HCOLOR— 2 
17210 HPLOT X2,Z2 

17220 IF Y2 < 50 OR Y2 > 191 THEN 17240 
17230 HPLOT X2,Y2 
17240 RETURN 
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APPENDIX B 


Sample Flutter Problem for Computer Program 

Due to the fact that errors can occur in entering the program in the 
computer, a sample case is given which can be used to determine if the program 
is working correctly. The sample problem appears in reference 4 pp. 219-224. 
The listed values for input are, 

u = 76 
a = -0.15 
x a = .25 

r 2 = .388 
a 

b = 5 inches 
aj a = 64.1 radians/sec 
= 55.9 radians/sec 
9h = 9ct = 0 


As most of these parameters are not in the units required as input for the 
program, the units are converted as follows: 


b = .4167 ft 

fy, = oj^/2ti = 8.9 Hz 

f 0 = co„/2ir = 10.2 Hz 

EA% = (1/2 + a/2)xl00 = 42.5% 

CG% = EA% + (x a /2 x 100) = 55% 


Assuming sea level air density 


p = .00237 slugs/ft3 
m = inrpb^ = .098 slugs/ft-span 


Therefore, 


I = r 2 mb 2 = .0066 slug-ft 2 /ft-span 

Furthermore, incompressible flow and infinite span are assumed. Therefore, 


M = 0 

AR = 100,000,000 (Essentially infinite) 

The input parameters appear in Table I. 

The full solution including the intermediate values is presented herein 
beginning on page 24. The solution is also plotted in Figure 6. The flutter 
solution is found by tracking one root of the quadratic to the crossing of the 
axis. The two successive k-values which have damping values on opposite sides 
of the axis are then averaged to determine the next k-value to be used. The 
solution process continues until the velocity difference for two successive 
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k-values is less than one foot per second. In other words, after first 
calculating that the flutter solution is between k-values of 0.3 and 0.2, the 
program then iterates to find that the flutter solution Is at a k-value of 
0.28125, The flutter velocity Is 89.3 ft/sec. The difference of this soluiton 
from the one found in reference 4 Is less than one percent. The reason for 
this difference may he due to the aerodynamic approximation of the program. 
The aerodynamics utilized by reference 4 are the actual values calculated from 
the Bessel functions whereas the algorithm used in the program is an 
approximation of these functions. 

It is also Interesting to note that the calculations at k-values of 0.01 
and 0.001 indicate a second Instability at about 173 ft/sec where both the 
frequency and damping of one of the modes goes to zero. This Is the static 
divergence speed. The divergence phenomena is discussed In textbooks such as 
references 3, 4, and 6, 
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MACH NUMBER s 0 AIR DENSITY a 2.37E-03 SLUGS/FT~3 
K-VALUE : 10 

FREQUENCIES (HZ) r 12.4688466 7.90074515 

VELOCITIES (F/S) s 3.26459755 2.06857571 

G-VALUE : -7. 63529253E-05 -1 . B7306349E-03 

K-VALUE s 6 

FREQUENCIES (HZ) s 12.4655681 . 7.90116384 

VELOCITIES (F/S) s 5.43956527 3.44780889 

G-VALUE a — 1 . 41403715E— 04 -3. 12198597E-03 

K-VALUE s 4 

FREQUENCIES (HZ) : 12.458B868 7.90205072 

VELOCITIES (F/S) a 8.15497465 5.17229385 

G-VALUE s -2. 52936026E— 04 -4. 683B3308E-03 

K-VALUE s 3 

FREQUENCIES (HZ) s 12.4494208 7.90331766 

VELOCITIES (F/S) s 10.8650383 6.8974975 

G-VALUE : -4.1 248 1 44BE-04 -6. 24684599E-03 

K-VALUE s 2 

FREQUENCIES (HZ) s 12. 4221688 7.90697101 

VELOCITIES (F/S) s 16.2618817 10.3510289 

G-VALUE s -9 . 33 1 0006 1 E-04 -9. 37824872E-03 

K-VALUE s 1.5 

FREQUENCIES (HZ) : 12.3835824 7.9121559 

VELOCITIES (F/S) s 21.6151576 13.8104219 

G-VALUE : -1 . 8052455E-03 -.0125208874 
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K-VALUE s 1.2 


FREQUENCIES (HZ) i 12.333302 
VELOCITIES (F/S) i 26.9092435 
G-VALUE I -3. 1 0577846E— 03 

K-VALUE s 1 

FREQUENCIES (HZ) i 12.2709408 
VELOCITIES (F/S) s 32.127818 
G-VALUE s -4.881 65B99E-03 


K-VALUE s .8 

FREQUENCIES (HZ) s 12.1539559 
VELOCITIES (F/S) s 39.7769097 
G-VALUE s —8. 45986823E— 03 

K-VALUE S .66 

FREQUENCIES (HZ) s 11.9981948 
VELOCITIES (F/S) : 47.5965354 
G-VALUE : -.0133534084 

K-VALUE : . 6 

FREQUENCIES (HZ) s 11.8942858 
VELOCITIES (F/S) : 51.902764 
G-VALUE s -.0165571556 

K-VALUE : .56 

FREQUENCIES (HZ) s 11.8050037 
VELOCITIES (F/S) s 55.1926781 
G-VALUE : -.0192229777 


7.91894171 

17.2778328 

-.0156811715 

7.9274217 

20.7556018 

-.0188676793 

7.94358842 

25.9974119 

-.0237219461 

7 . 9658^205 
31.600175 
.0290210464 

7.98119105 

34.8273015 

.0321504883 

7.99481965 

37.3786845 

.0346779895 
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K-VALUE 


5 


FREQUENCIES (HZ) s 11.6273907 S. 02324407 

VELOCITIES (F/S) s 60.8857458 42.0129684 

G-VALUE s 0242034231 0394240991 

K-VALUE s .4 

FREQUENCIES <HZ> : 11.1289104 8.11506943 

VELOCITIES (F/S> : 72.8443754 53.1172543 

G-VALUE s -.0351844704 -.052067462 

K-VALUE : .3 

FREQUENCIES (HZ) s 9.99785211 8.44082058 

VELOCITIES (F/S) : 87.2546984 73.6659481 

G-VALUE : -.0347203235 -.0898817964 

K-VALUE s . 2 

FREQUENCIES (HZ) s 8.81288397 8.56418574 

VELOCITIES (F/S) : 115.36961 112.113898 

G-VALUE s .243699161 -.437955103 

K-VALUE s .16 

FREQUENCIES (HZ) : 8.4513254B 8.12400021 

VELOCITIES (F/S) s 138.295552 132.939277 

G-VALUE : .332460844 -.561335102 

K-VALUE s .12 

FREQUENCIES (HZ) s 7.85576763 7.2834035 

VELOCITIES (F/S) : 171.399974 158.911927 

G-VALUE : .382571844 -.632265805 
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K-VALUE 8 .1 


FREQUENCIES (HZ) : 7.42124838 6.5961167 

VELOC ITIES (F/S) ; 194. 30337 1 172. 699746 

G-VALUE J . 366 175506 -.61 0284063 

K-VALUE : .08 

FREQUENCIES (HZ) : 6.93017177 5.60964015 

VELOCITIES (F/S) x 226.807484 183.58973 

G-VALUE s .285409318 -.506082019 

K-VALUE : . 06 

FREQUENCIES (HZ) x 6.69811934 4.23294678 

VELOCITIES (F/S) s 292.283971 184.711921 

G-VALUE s .1471 28775 - . 334656658 

K-VALUE : .04 

FREQUENCIES (HZ) s 6. 76022733 2. 75350606 

VELOCITIES (F/S) s 442.49124 180. 230967 

G-VALUE s . 066832 1 095 -.21 9908692 

K-VALUE : .025 

FREQUENCIES (HZ) r 6.80672942 1.68681855 

VELOCITIES (F/S) s 712.856061 176.657357 

G-VALUE s .0351797178 -. 147910484 

K-VALUE s .01 

FREQUENCIES (HZ) s 6.83071834 .66366487 

VELOCITIES (F/S) s 1788.42094 173.760956 

G-VALUE : .01282117 -.0636466531 
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K- VALUE s IE-03 


FREQUENCIES (HZ) s 6.57653004 .066106868B 

VELOCITIES (F/S) s 17218.6928 173.081223 

G-VALUE s 1 . 64635228E— 03 -6 . 47606B03E-03 

K-VALUE s . 25 

FREQUENCIES s 9.13877632 8.80154616 

VELOCITIES • 95.7086978 92.1769493 

G-VALUES s .102772118 -.256309331 

K-VALUE s . 275 

FREQUENCIES s 9.44658006 8.70630525 

VELOCITIES s 89.9384262 82.8904628 

G-VALUES s 5. 4208891 5E— 03 -.142427175 

K-VALUE s . 2875 

FREQUENCIES s 9.7390505 8.55501294 

VELOCITIES : 88.691528 77.9087416 

G-VALUES s -.0240007964 -.106496951 

K-VALUE s .28125 

FREQUENCIES : 9.59222603 8.62878413 

VELOCITIES s 89.2956372 80.3267954 

G-VALUES s -.012886655 -.120744428 


THE AIRCRAFT WILL FLUTTER NEAR THE FOLLOWING VALUES 
K-VALUE : .28125 

FREQUENCIES : 9.59222603 8.62878413 

VELOCITIES : 89.2956372 80.3267954 

Note: When the program is executed, 

g-values : -.012886655 -.120744428 a plot of the results similar 

to the one shown in Figure 6 
will be output following the 
listing. 
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APPENDIX C 


A matched point solution is presented in this section. The solution is 
calculated for four different Mach numbers and the solution Is presented In 
terms of dynamic pressure, altitude and equivalent airspeed. 

The wing to be analyzed is taken from reference 3, page 203. An aspect 
ratio of eight was assumed. The Input appears in table II. Air density 
option two was chosen and calculations were made for six value of air 
density. The interval between values was 0.0004 and the starting value was 
0.0004. The speed of sound at each air density was obtained from a standard 
atmosphere chart and then multiplied by the appropriate Mach number to obtain 
the corresponding velocity. 

The results for each Mach number are in Table III* 

Results at M = .4 . When the flutter velocity is plotted versus air density as 
shown in Figure 7, it is seen that the aircraft does not flutter at M = .4 at 
any altitude above sea level. Therefore, at M = .4 the aircraft is safe from 
flutter. 

Results at M = .5 . The velocity is plotted versus air density in Figure 8 the 
curve intersects the Mach 0.5 curve at an air density value of .0021 
slugs/ft 3 and a velocity of 550 ft/sec. 

Therefore, 

q = l/2pV 2 = 317.6 lb/ft 2 


The altitude H, In feet, can be calculated for altitudes less than about 
40,000 feet using the following equation which was derived from equations 
obtained from reference 5. 


H = 14S45O.[l-{- TJ0 § 37 .) 


.23496, 


Thus for p = .0021 the altitude is H = 4075 ft. 
The equivalent airspeed is: 


V e * V-3TO23T = 517 - 7 ft/sec 


Therefore, the aircraft is safe from flutter at M = .5 for altitudes greater 
than 4075 ft. If the plane was at 4075 ft., the aircraft would flutter when 
the equivalent airspeed is 517 ft/sec. For a more detailed description of the 
relevance of equivalent airspeed see reference 8. 
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Results at M = .6 


Using the above procedure, a matched point is found by using Figure 9. The 
plot indicates the curves crossing at an air density of 0.0014 slugs/ft 3 and 
a velocity of 630 ft/sec. The flutter conditions are therefore: 


q = 277.8 lb/ft2 
H = 16922 ft 
V e = 484.2 ft/sec 


Results at M = .8 

Using Figure 10, the matched point occurs at 

p = .00065 slugs/ft 3 
V = 780 ft/sec 


Therefore, 


q = 197.8 lbs/ft 2 
H = 38124 ft 
V e = 408.5 ft/sec 


When plotted as dynamic pressure versus Mach number as shown in Figure 11, the 
flutter boundary is described. However, when plotted as equivalent airspeed 
versus Mach number as in Figure 12, the flutter boundary takes on added 
meaning when the altitude lines are drawn. This plot makes it easy to 
visualize where the aircraft is safe in terms of altitude and equivalent air 
speed. The abbreviation "K" means thousand and is utilized in terms of feet 
( i . e. , 10K means 10,000 feet altitude). 
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Table I.- Input for problem in Appendix B 


INPUT PARAMETERS : 

MASS PER FOOT SPAN : .098 SLUGS 

CG LOCATION s 55 7. CHORD 

PITCH INERTIA s 6.6E-03 SLUG-FT^ 

ELASTIC AXIS LOCATION s 42.5 7. CHORD 
SEMI -CHORD : .4167 FEET 

NATURAL BENDING FREQUENCY : 8.9 HZ 
NATURAL TORSION FREQUENCY 8 10.2 HZ 
BENDING DAMPING : O 7. CRITICAL DAMPING 
TORSION DAMPING s 0 7.CRITICAL DAMPING 
FULL SPAN ASPECT RATIO : 100000000 

MAXIMUM VELOCITY TO BE PLOTTED s 200 FEET/SEC 
MAXIMUM FREQUENCY TO BE PLOTTED : 15 HZ 
MAXIMUM G-VALUE TO BE PLOTTED : .05 
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Table II.- Input for problem in Appendix C. 


INPUT PARAMETERS : 

MASS PER FOOT SPAN : .6516 SLUGS 
CG LOCATION s 46 V. CHORD 
PITCH INERTIA : 3.375 SLUG-FT^ 

ELASTIC AXIS LOCATION : 35 % CHORD 
SEMI -CHORD : 3.125 FEET 
NATURAL BENDING FREQUENCY s 9.9 HZ 
NATURAL TORSION FREQUENCY : 16.02 HZ 
BENDING DAMPING : 0 7. CRITICAL DAMPING 
TORSION DAMPING : O ‘/.CRITICAL DAMPING 
FULL SPAN ASPECT RATIO : 8 

MAXIMUM VELOCITY TO BE PLOTTED : 1200 FEET /SEC 
MAXIMUM FREQUENCY TO BE PLOTTED : 20 HZ 
MAXIMUM G-VALUE TO BE PLOTTED : .05 
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Table III.- Results of Matched Point Problem in Appendix C 


Mach Number 

Air Density 

Flutter Velocity 

Mach Number 

M 

P, o 

V, 

Velocity, MV , ft/sec 


slugs/ft 0 

ft/sec 

a 

.4 

.0004 

1192 

387 


.0008 

862 

393 


.0012 

717 

412 


.0016 

630 

426 


.0020 

572 

438 


.0024 

529 

447 

.5 

.0004 

1161 

484 


.0008 

840 

492 


.0012 

699 

515 


.0016 

614 

533 


.0020 

559 

547 


.0024 

516 

559 

.6 

.0004 

1118 

581 


.0008 

810 

590 


.0012 

674 

618 


.0016 

594 

640 


.0020 

539 

656 


.0024 

500 

670 

.8 

.0004 

977 

774 


.0008 

709 

786 


.0012 

593 

824 


.0016 

525 

853 


.0020 

480 

875 


.0024 

448 

894 
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Damping Frequency 




b b H 
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Figure 4.- Location of center of gravity of a wing. 



Figure 5.- Location of elastic axis of a wing. 
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Mach number: 0 Air density: 2.37E-03 slugs/ft 
The aircraft will flutter near the following values 
K-value: .28125 

Frequencies: 9.59222603 8.62878413 

Velocities : 89 . 2956372 80.3267954 

'G- values : -.012886655 -.120744428 



Figure 6.- Solution of sample problem in Appendix B. 
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Figure 9.- Matched point for Mach .6. 
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Figure 12.- Equivalent airspeed flutter boundary. 
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